/**
 * mLibMat44.js
 * Created: Nov 2013
 * Source link: http://code.google.com/p/code-depository/source/browse/
 *
 * By Suresh Nalluri,
 * License: The MIT License (MIT)
 *          See LICENSE details
 *
 *
 *	Todo:
 *
 *
 * Requires:
 *	mLibIncludes.js
 *  mLibMat33.js
 *	mLibFactory.js
 *	mLibVec3.js
 *	mLibVec4.js
 *	mLibQuaternion.js,
 *
 * 4x4 Matrix Class
 * A matrix object is composed of a one dimensional array of 9 elements.  The matrix is sorted in
 * column major order, which imples that element [i,j] is the (i+4*j)th entry in the array.
 *
 * observe that the array rep of the matrix stores the matrix column with following index:
 *
 * 0 4  8 12
 * 1 5  9 13
 * 2 6 10 14
 * 3 7 11 15
 *
 * eg. col 1 of  the matrix is the transpose of: ( A[0], A[1], A[2], A[3] )
 * where A is the array storing the matrix.
 *
 * note this constructor pattern allows us to create objects in 2 ways:
 * 1. var o = mLib.Mat44( 0 or more params );      //OR
 * 2. var o = new mLib.Mat44( 0 or more params );
 *
 * NOTE: var o = [new] mLIB.Mat44(<no params>); // will return default identity matrix *
 *
 * added benefit of this pattern is that we get a copy constructor by default!!
 * eg. var newMtx = mLIB.Mat44(existingMtx);
 * so no need to have an explicit copy method.
 *
 * Usage:
 *		var l = [new] mLIB.Mat44( mat | quat |  <noargs> ]);
 * 
 *-------------------------------------------------------
 **/
if (typeof mLIB == 'undefined') function mLIB() {};

( function(module) {

    var Vec3 = module.Vec3,
        Vec4 = module.Vec4,
        
        eqZero = module.Includes.isZero,
        sqrt = module.Includes.sqrt,
        sin = module.Includes.sin,
        cos = module.Includes.cos,
        acos = module.Includes.acos,
        atan2 = module.Includes.atan2,
        PI = module.Includes.PI,
        ZERO = module.Includes.ZERO,
        axisIdx = ['x', 'y', 'z'];

	module.Mat44 = Factory.obj( {
        init: function() {
            var A = this.mat = new Array(16);
            
            if (arguments[0].isClass(module.Mat44) ) {
                // param is a matrix - copy constructor
                var src = arguments[0].mat;
                
                A[0] = src[0]; A[4] = src[4]; A[8]  = src[8];  A[12] = src[12];
                A[1] = src[1]; A[5] = src[5]; A[9]  = src[9];  A[13] = src[13];
                A[2] = src[2]; A[6] = src[6]; A[10] = src[10]; A[14] = src[14];
                A[3] = src[3]; A[7] = src[7]; A[11] = src[11]; A[15] = src[15];
                
            } else if (arguments[0].isClass(module.Quat) {
                // param is a quaternion - axis-angle construction
                this.rotationByQuaternion(arguments[0]);
                
            } esle if (arguments[0].isClass(module.Mat33) {
                // param is a 3x3 matrix.. use this.rotationByMat33(src)
                this.rotationByMat33(arguments[0]);
                
            } else {
                // identity matrix
                A[0] = 1.0; A[4] = 0.0; A[8]  = 0.0; A[12] = 0.0;
                A[1] = 0.0; A[5] = 1.0; A[9]  = 0.0; A[13] = 0.0;
                A[2] = 0.0; A[6] = 0.0; A[10] = 1.0; A[14] = 0.0;
                A[3] = 0.0; A[7] = 0.0; A[11] = 0.0; A[15] = 1.0;
                  
            }
        }, //init()
        
        // synonym for the '=' operand, eg. mtx.equals(other) - copies other into mtx.
        equals: function(m) {
            if ( this === m ) return;
            
            var A = this.mat,
                B = M.mat;
            
            A[0] = B[0]; A[3] = B[3]; A[6] = B[6];
            A[1] = B[1]; A[4] = B[4]; A[7] = B[7];
            A[2] = B[2]; A[5] = B[5]; A[8] = B[8];
            A[9] = B[9];
            A[10] = B[10];
            A[11] = B[11];
            A[12] = B[12];
            A[13] = B[13];
            A[14] = B[14];
            A[15 = B[15];
            
            return this;
        },

        CopyOf: function(m) {
            // ie. clone
            return module.Mat44( m );
        },
        
        // clean self.. set near-zero elements to zero.
        clean: function() {
            var A = this.mat;
            
            for (var i = 0; i < 15; i++) {
                eqZero(A[i]) && (A[i] = 0.0);
            }
        },
        
        Clean: function() {
            return module.Mat44(this).zero();
        },
        
        // set self to identity matrix
        identity: function() {
            var A = this.mat;
            
            A[0] = 1.0; A[4] = 0.0; A[8]  = 0.0; A[12] = 0.0;
            A[1] = 0.0; A[5] = 1.0; A[9]  = 0.0; A[13] = 0.0;
            A[2] = 0.0; A[6] = 0.0; A[10] = 1.0; A[14] = 0.0;
            A[3] = 0.0; A[7] = 0.0; A[11] = 0.0; A[15] = 1.0;
            return this;
        },
       

        // comparison operations
        
        isEqual: function(m) {
            if ( this === m )return true;  // since it's the same matrix object!!
            
            var A = this.mat,
                B = m.mat;
                
            for (var i = 0; i < 16; i++) {
                if ( !eqZero(A[i] - B[i]) return false;
            }
            return true;
        },
        
        isNotEqual: function(m) {
            var A = this.mat,
                B = m.mat;
                
            for (var i = 0; i < 16; i++) {
                if ( !eqZero(A[i] - B[i]) return true;
            }
            return false;
       },
       
        isZero: function() {
            var A = this.mat;                

            for (var i = 0; i < 16; i++) {
                if ( !eqZero(A[i]) return false;
            }
            return true;
        },
        
        isIdentity: function() {
            var A = this.mat;
            
            return (
                eqZero(1.0-A[0]) && eqZero(A[4]) && eqZero(A[8])  && eqZero(A[12]) &&  
                eqZero(A[1]) && eqZero(1.0-A[5]) && eqZero(A[9])  && eqZero(A[13]) &&  
                eqZero(A[2]) && eqZero(A[6]) && eqZero(1.0-A[10]) && eqZero(A[14]) &&  
                eqZero(A[3]) && eqZero(A[7]) && eqZero(A[11]) && eqZero(1.0-A[15])
            );
        
        },
        
        // transpose self
        transpose: function() {
            var tmp,
                A = this.mat;
            
            tmp = A[4];
            A[4] = A[1];
            A[1] = tmp;
            
            tmp = A[8];
            A[8] = A[2];
            A[2] = tmp;
            
            tmp = A[9];
            A[9] = A[6];
            A[6] = tmp;
            
            tmp = A[12];
            A[12] = A[3];
            A[3] = tmp;
            
            tmp = A[13];
            A[13] = A[7];
            A[7] = tmp;
            
            tmp = A[14];
            A[14] = A[11];
            A[11] = tmp;

            return this;
        },

        // return the transpose of self in a new Mat44
        Transpose: function() {
            return new module.Mat44(this).transpose();
        }
        
        // Calculate the inverse of a standard 4x4 affine matrix.
        // An affine matrix is one where the bottom row has the form
        // [0 0 ... 0 1]
        // We use Cramer's rule to calculate the inverse as it is much faster than
        // using gaussian elimination for a 4x4.
        // Cramer's rule states that inv(A) = (1/det(A))adj(A)
        // where inv(A) means inverse of A and adj(A) means the adjoint matrix of A
        // recall that the adjoint is the TRANSPOSE of the matrix composed of the
        // cofactors C[ij] of each element of A.  where C[ij] = (-1)^(i+j)*det(sub(A[ij])
        // sub(A[ij]) is the sub-matrix of A, formed by removing the the row i and
        // column j from A.  Therefore adj(A) = transpose(C[ij])
        //
        affineInverse: function() {
            // set self as the inverse
            // so our 4x4 looks as follows:
            //                      note the signs for cofactors
            //      a b c d             + - + -
            //      e f g h             - + - +
            //      i j k l             + - + -
            //      0 0 0 1             - + - +
            //
            //  since bottom row has one non zero colum, we use this row ro compute det(A)
            //  det(A) = det(upper left 3x3 matrix)
            // so lets get cofactors for this 3x3 using row 0            
            var A = this.mat,
                cf0 = A[5]*A[10] - A[6]*A[9],       // fk-gj
                cf4 = A[2]*A[9] - A[1]*A[10],       // -(ek-gi) = gi - ek
                cf8 = A[1]*A[6] - A[2]+A[5],        //ej-fi
                det = A[0]*cf0 + A[4]*cf4 + A[8]*cf8; // a*cf0 + b*cf4 + c *cf8
            
            if (eqZero(det)) return null;
            
            // let's create the rest of the cofactors of A now
            // the cofactor matrix is as follows:
            //      cf0 cf4 cf8  0       and the adjoint is:    cf0 cf1 cf2  0
            //      cf1 cf5 cf9  0                              cf4 cf5 cf6  0
            //      cf2 cf6 cf10 0                              cf8 cf9 cf10 0
            //      0   0   0    det                            0   0   0    det
            //
            // where:
            //  cf1 = -(bk - cj) = cj - bk
            //  cf5 =  ak - ci
            //  cf9 = -(aj - bi) = bi - aj
            //  cf2 = bg - cf
            //  cf6 = -(ag-ce) = ce - ag
            //  cf10 = af - be
            
            var recipDet = 1.0 / det;                
            
            // set A's upper left 3x3 to the transpose of the cofactors / det. ie the adjoint/det
            A[0] = recipDet*cf0;
            A[1] = recipDet*cf4;
            A[2] = recipDet*cf8;
            A[4] = recipDet*(A[6]*A[8] - A[4]*A[10]);   // cf1
            A[5] = recipDet*(A[0]*A[10] - A[2]*A[8]);   // cf5
            A[6] = recipDet*(A[2]*A[4] - A[0]*A[6]);    // cf9
            A[8] = recipDet*(A[4]*A[9] - A[5]*A[8]);    // cf2
            A[9] = recipDet*(A[1]*A[8] - A[0]*A[9]);    // cf6
            A[10] = recipDet*(A[0]*A[5] - A[1]*A[4]);   // cf10
                
            // Note cf3 = cf7 = cf11 = 0.0 = cf12 = cf13 = cf14 = 0.0 so no need to worry about them
            // c15 = det; ie the determinant we calculated earlier
                
            // in block matrix form, our Affine matrix A looks as follows:
            //      M   b
            //      0   1
            // 
            //  therefore if inv(A) =  X c
            //                         0 1     then
            //
            // Ainv(A) = MX Mc+b    = I 0
            //           0  1         0 1   which means that
            //
            // MX = I => X = inv(M)  and Mc + b = 0 => Mc = -b => c = -inv(M)b
            // we have inv(M) above so final step is to evaluate c
            //
            var t12 = A[12], t13 = A[13], t14 = A[14];  // the b vector
            A[12] = -A[0]*t12 - A[4]*t13 - A[8]*t14;
            A[13] = -A[1]*t12 - A[5]*t13 - A[9]*t14;
            A[14] = -A[2]*t12 - A[6]*t13 - A[10]*t14;
            
            // note the bottom row of A should by default be [0 0 0 1], since A was an
            // affine matrix to start with!
            
            return this;
        },
        
        AffineInverse: function(m) {
            // return inverse in new matrix object
            return Mat44(m).affineInverse();
        },
        
        multiplyByScalar: function(f) {
            var A = this.mat;
            A[0] *= f; A[4] *= f; A[8]  *= f; A[12] *= f;
            A[1] *= f; A[5] *= f; A[9]  *= f; A[13] *= f;
            A[2] *= f; A[6] *= f; A[10] *= f; A[14] *= f;
            A[3] *= f; A[7] *= f; A[11] *= f; A[15] *= f;
            
            return this;
        },
        
        MultiplyByScalar: function(f) {
            return new Mat44(this).multiplyByScalar(f);
        },
        
        plus: function(m) {
            var A = this.mat,
                B = m.mat;
                
            for (var i =0; i< 16; i++) {
                A[i] += B[i];
            }
            return this;
        },
        
        Plus: function(m) {
            return module.Mat44(this).plus(m);
        },
        
        minus: function(m) {
            var A = this.mat,
                B = m.mat;
                
            for (var i =0; i< 16; i++) {
                A[i] -= B[i];
            }
            return this;
        },

        Minus: function(m) {
            return module.Mat44(this).minus(m);
        },

        negation: function() {
            var A = this.mat;
            
            for (var i =0; i< 16; i++) {
                A[i] = -A[i];
            }
            return this;
        
        },
        
        Negation: function() {
            return module.Mat44(this).negation();
        },
        
        multiply: function(m) {
            var A = this.toArray(),  // cheap than creating a copy of this?
                B = m.mat,
                C = this.mat;
            
            C[0]  = A[0]*B[0]  + A[4]*B[1]  + A[8]*B[2]  + A[12]*B[3];
            C[4]  = A[0]*B[4]  + A[4]*B[5]  + A[8]*B[6]  + A[12]*B[7];
            C[8]  = A[0]*B[8]  + A[4]*B[9]  + A[8]*B[10] + A[12]*B[11];
            C[12] = A[0]*B[12] + A[4]*B[13] + A[8]*B[14] + A[12]*B[15];

            C[1]  = A[1]*B[0]  + A[5]*B[1]  + A[9]*B[2]  + A[13]*B[3];
            C[5]  = A[1]*B[4]  + A[5]*B[5]  + A[9]*B[6]  + A[13]*B[7];
            C[9]  = A[1]*B[8]  + A[5]*B[9]  + A[9]*B[10] + A[13]*B[11];
            C[13] = A[1]*B[12] + A[5]*B[13] + A[9]*B[14] + A[13]*B[15];
            
            C[2]  = A[2]*B[0]  + A[6]*B[1]  + A[10]*B[2]  + A[14]*B[3];
            C[6]  = A[2]*B[4]  + A[6]*B[5]  + A[10]*B[6]  + A[14]*B[7];
            C[10] = A[2]*B[8]  + A[6]*B[9]  + A[10]*B[10] + A[14]*B[11];
            C[14] = A[2]*B[12] + A[6]*B[13] + A[10]*B[14] + A[14]*B[15];
            
            C[3]  = A[3]*B[0]  + A[7]*B[1]  + A[11]*B[2]  + A[15]*B[3];
            C[7]  = A[3]*B[4]  + A[7]*B[5]  + A[11]*B[6]  + A[15]*B[7];
            C[11] = A[3]*B[8]  + A[7]*B[9]  + A[11]*B[10] + A[15]*B[11];
            C[15] = A[3]*B[12] + A[7]*B[13] + A[11]*B[14] + A[15]*B[15];
            
            return this;
        },
        
        Multiply: function(m) {
            return new Mat44(this).multiply(m);
        },
        
        // Matrix - Column vector multiply ie c = [A]b
        multiplyColumn: function(colVec) {
            var A = this.mat,
                b = colVec.toArray();       // this may not be cheap !!
            
            return module.Vec4( 
                A[0]*b[0] + A[4]*b[1] + A[8]*b[2]  + A[12]*b[3],
                A[1]*b[0] + A[5]*b[1] + A[9]*b[2]  + A[13]*b[3],
                A[2]*b[0] + A[6]*b[1] + A[10]*b[2] + A[14]*b[3],
                A[3]*b[0] + A[7]*b[1] + A[11]*b[2] + A[15]*b[3]
            );
        },    
        
        // Row vector - Matrix multiply  ie. c = b[A]
        multiplyRow: function(rowVec) {
            var A = this.mat,
                b = rowVec.toArray();
            
            return module.Vec4(
                b[0]*A[0]  + b[1]*A[1]  + b[2]*A[2]  + b[3]*A[3],
                b[0]*A[4]  + b[1]*A[5]  + b[2]*A[6]  + b[3]*A[7],
                b[0]*A[8]  + b[1]*A[9]  + b[2]*A[10] + b[3]*A[11],
                b[0]*A[12] + b[1]*A[13] + b[2]*A[14] + b[3]*A[15]
            );
        },
        
        // Set translation matrix, where translation is given as a 3D vector t
        translation: function(t) {
            var A = this.mat;
            
            A[0] = A[5] = A[10] = A[15] = 1.0;
            A[1]=A[2]=A[3]=A[4]=A[6]=A[7]=A[8]=A[9]=A[11]=0.0;
            A[12] = t.x;
            A[13] = t.y;
            A[14] = t.z;

            return this;
        },
        
        addTranslation: function(xlation) {
            // this is a handy compromise method to add translation quickly
            var A = this.mat;
            
            A[12] = xlation.x;
            A[13] = xlation.y;
            A[14] = xlation.z;
            
            return this;
        },
        

        rotationByQuaternion: function(Q) {
        
            var xs = Q.x+Q.x,
                ys = Q.y+Q.y,
                zs = Q.z+Q.z,
                wx = Q.w*xs,
                wy = Q.w*ys,
                wz = Q.w*zs,
                xx = Q.x*xs,
                xy = Q.x*ys,
                xz = Q.x*zs,
                yy = Q.y*ys,
                yz = Q.y*zs,
                zz = Q.z*zs,
                A = this.mat;

            A[0] = 1.0 - (yy + zz);
            A[4] = xy - wz;
            A[8] = xz + wy;
            A[12] = 0.0;
            
            A[1] = xy + wz;
            A[5] = 1.0 - (xx + zz);
            A[9] = yz - wx;
            A[13] = 0.0;
            
            A[2] = xz - wy;
            A[6] = yz + wx;
            A[10] = 1.0 - (xx + yy);
            A[14] = 0.0;

            A[3] = 0.0;
            A[7] = 0.0;
            A[11] = 0.0;
            A[15] = 1.0;

            return *this;
        
        },
        
        // rotation based on a 3x3 rotation matrix
        rotationByMat33: function(m3) {
            var A = this.mat, M = m3.mat;
            
            A[0]  = B[0];
            A[1]  = B[1];
            A[2]  = B[2];
            A[3]  = 0.0;
            A[4]  = B[3];
            A[5]  = B[4];
            A[6]  = B[5];
            A[7]  = 0.0;
            A[8]  = B[6];
            A[9]  = B[7];
            A[10] = B[8];
            A[11] = A[12] = A[13] = A[14] = 0.0;
            A[15] = 1.0;
            
            return this;
            
        },
                
        // rotate around a given axis, where axis denotes a vector around which we rotate
        //
        rotationByAxisAngle: function(axis, angle) {
            var c = cos(angle), s = sin(angle), p = 1.0 - c,
                n = axis.Normalise(),       // normalised vector in the axis direction
                px = p * n.x, py = p * n.y, pz = p * n.z,
                sx = s * n.x, sy = s * n.y, sz = s * n.z,
                pxy = px * n.y, pyz = py * n.z, pxz = px * n.z,
                A = this.mat;
                
            A[0] = px*n.x + c; A[4] = pxy - sz;   A[8] = pxz + sy;
            A[1] = pxy + sz;   A[5] = py*n.y + c; A[9] = pyz - sx;
            A[2] = pxz - sy;   A[6] = pyz + sx;   A[10] = pz*n.z + c;
            A[3] = A[7] = A[11] = A[12] = A[13] = A[14] = 0.0;
            A[15] = 1.0;

            return this;
        },
        
        // set self as rotation matrix, based on angle in radians around the x-axis
        rotationX: function(theta) {
            var sinTheta = sin(theta),
                A = this.mat;
            
            // +-                -+
            // | 1.0 0.0  0.0 0.0 |
            // | 0.0 cos -sin 0.0 |
            // | 0.0 sin  cos 0.0 |
            // | 0.0 0.0  0.0 1.0 |
            // +-                -+
            
            A[0] = A[15] = 1.0;
            A[1] = A[2] = A[3] = A[4] = 0.0;
            A[5] = A[10] = cos(theta);
            A[6] = sinTheta;
            A[7] = A[8] = 0.0;
            A[9] = -sinTheta;
            A[11] = A[12] = A[13] = A[14] = 0.0;
            
            return this;
        },
        
         // set self as rotation matrix, based on angle in radians around the y-axis
        rotationY: function(theta) {
            var sinTheta = sin(theta),
                A = this.mat;

            // +-                -+
            // | cos  0.0 sin 0.0 |
            // | 0.0  1.0 0.0 0.0 |
            // | -sin 0.0 cos 0.0 |
            // | 0.0  0.0 0.0 1.0 |
            // +-                -+

            A[0] = A[10] = cos(theta);
            A[2] = -sinTheta;
            A[5] = A[15] = 1.0;
            A[8] = sinTheta;
            A[1] = A[3] = A[4] = A[6] = A[7] = A[9] = A[11] = A[12] = A[13] = A[14] = 0.0;
            A[15] = 1.0;
            
            return this;
        },
         // set self as rotation matrix, based on angle in rdians around the z-axis
        rotationZ: function(radians) {
            var sinTheta = sin(theta),
                A = this.mat;

            // +-                -+
            // | cos -sin 0.0 0.0 |
            // | sin  cos 0.0 0.0 |
            // | 0.0  0.0 1.0 0.0 |
            // | 0.0  0.0 0.0 1.0 |
            // +-                -+
            
            A[0] = A[5] = cos(theta);
            A[1] = sinTheta; A[4] = -sinTheta;
            A[10] = A[15] = 1.0;
            A[2] = A[3] = A[6] = A[7] = A[8] = A[9] = A[11] = A[12] = A[12] = A[14] = 0.0;
            
            return this;
        },

        // set the matrix to a rotation matrix - given the rotation around each axis individually
        // ie rotation based on euler angle
        //
        rotationByEuler: function(xAngle, yAngle, zAngle) {
            // if the rotations by axis are represented by the matrices: Rx, Ry, Rz, then
            // the rotation = RxRyRz
            var cx = cos(xAngle), sx = sin(xAngle),
                cy = cos(yAngle), sy = sin(yAngle),
                cz = cos(zAngle), sz = sin(zAngle),
                A = this.Mat;
            
            A[0] =  cy*cz;            A[4] = -cy*sz;            A[8] = sy;     A[12] = 0.0;
            A[1] =  sx*sy*cz + cx*sz; A[5] = -sx*sy*sz + cx*cz; A[9] = -sx*cy; A[13] = 0.0;
            A[2] = -cx*sy*cz + sx*sz; A[6] =  cx*sy*sz + sx*cz; A[10] = cx*cy; A[14] = 0.0;
            A[3] = A[7] = A[11] = 0.0; A[15] = 1.0;
            
            return this;
        },

        // get one set of possible z-y-x fixed angles that will generate this matrix
        // Assumes that this is the fixed angles (z-x-y) rotation matrix:
        //  +-                                           -+
        //  |  cy*cz             -cy*sz           sy    0 |   where cy = Cos(yTheta)
        //  |  sx*sy*cz+cx*sz    -sx*sy*sz+cz*cz -sx*cy 0 |
        //  | -cx*sy*cz+sx*sz     cx*sy*sz+sx*cz  cx*cy 0 |
        //  |  0                  0               0     1 |
        //  +-                                           -+
        //
        // This is also the Euler angles (roll-pitch-yaw or x-y-z) matrix.
        //
        // note that what we return is one possible solution out of many due to
        // using inverse trig functions, which can result in multiple results
        //
        getFixedAngles: function(zyx) {
            var sx, sy, sz,
                cx, cy, cz,
                A = this.mat;
            
            sy = A[8];
            cy = sqrt(1.0 - sy*sy);
            if ( eqZero(cy) ) {
                // we have gimbal lock problem: x and z axes are aligned, which implies
                // we can't tell which is rotation around x and which is rotation around
                // y.  So let's assume that rotation around z is zero
                sz = 0.0;
                cz = 1.0;
                sx = A[6];
                cx = A[5];
                
            } else {
                // no gimbal lock
                var f = 1.0 / cy;
                sx = -A[9]*f;
                cx = A[10]*f;
                sz = -A[4]*f;
                cz = A[0]*f;
                
            }
            
            // return the angles in the given array param.
            // note use of atan2 as it returns angles within the range: [-PI, PI]
            zyx[0] = atan2(sz, cz); // the zTheta
            zyx[1] = atan2(sy, cy); // the yTheta
            zyx[2] = atan2(sz, cz); // the xTheta
        
            return this;
        },

        // Gets one possible axis-angle pair that will generate this matrix
        // Assumes that this is the rotation matrix around an arbitrary axis r:
        //  +-                                       -+
        //  | p*x^2+c       p*x*y-s*z   p*x*z+s*y   0 |   p = 1 - cos(theta)
        //  | p*x*y+s*z     p*y^2+c     p*y*z-s*x   0 |   c = cos(theta), s = sin(theta)
        //  | p*x*z-s*y     o*y*z+s*x   t*z^2+c     0 |   rAxis = (x, y, z)
        //  | 0             0           0           1 |
        //  +-                                       -+
        //  based on http://www.geometrictools.com/Documentation/RotationIssues.pdf
        //
        getAxisAngle: function(rAxis) {
            // we return the axis in the given Vec3 type parameter rAxis
            // and the angle theta in the return statement.
        
            var A = this.mat,
                trace = A[0] + A[5] + A[10],     // = 2Cos(theta) + 1
                theta = acos( (trace - 1.0) * 0.5 );
            
            if (eqZero(theta)) {
                // theta is zero, so we can return any old axis!
                // in our case we return the x Axis
                rAxis.x = module.Vec3.xAxis.x;
                rAxis.y = module.Vec3.xAxis.y;
                rAxis.z = module.Vec3.xAxis.z;
                return 0.0;
            }

            if (theta < PI-ZERO) {
                // our standard case, theta lies in range (0, PI)
                rAxis.x = A[6] - A[9];
                rAxis.y = A[8] - A[2];
                rAxis.z = A[1] - A[4];
                return theta;
            }
            
            // theta == PI
            var i = 0;
            
            if ( A[5] > A[0] )
                i = 1;
            if ( A[10] > A[i + 4*i] )
                i = 2;
            
            var j = (i+1)%3,
                k = (j+1)%3,
                s = sqrt( A[i + 4*i] - A[j + 4*j] - A[k + 4*k] + 1.0 ),
                f = 1.0 / s;
            
            rAxis[axisIdx[i]] = 0.5 * s;
            rAxis[axisIdx[j]] = (A[i + 4*j])*f;
            rAxis[axisIdx[k]] = (A[k + 4*i])*f;

            return theta;
        },
        
        // set self as scaling matrix based on a vector of scaling factors
        // ie. set the diagnonal elements of the matrix to the vector and the rest as 0's 
        scaling: function(sf) {
            var A = this.mat;
            
            A[0] = sf.x;
            A[5] = sf.y;
            A[10] = sf.z;
            A[1]=A[2]=A[3]=A[4]=A[6]=A[7]=A[8]=A[9]=A[11]=A[12]=A[13]=A[14] = 0.0;
            A[15] = 1.0f;
            return this;
        },

        // matrix - point multiplation
        Transform: function(point) {
            var A = this.mat,
                p = point.toArray();
                
            return Vec3 (
                A[0]*p[0] + A[4]*p[1] + A[8]*p[2]  + A[12],
                A[1]*p[0] + A[5]*p[1] + A[9]*p[2]  + A[13],
                A[2]*p[0] + A[6]*p[1] + A[10]*p[2] + A[14],
            );
        },
        
        TransformPoint: function(point) {   // same as transform method
            var A = this.mat,
                p = point.toArray();
                
            return Vec3(
                A[0]*p[0] + A[4]*p[1] + A[8]*p[2]  + A[12],
                A[1]*p[0] + A[5]*p[1] + A[9]*p[2]  + A[13],
                A[2]*p[0] + A[6]*p[1] + A[10]*p[2] + A[14],
            );
        },
        
        /**
         * Accessor methods.. setters and getters
         **/
        getElement: function(iRow, jCol) {
            return this.mat[iRow + 4*jCol];
        },
        
        $reader: function() {
            return ( function(iRow, jCol) {
                return this.mat[iRow + 4*jCol];
            }).bind(this);
        },
        
        // return the raw matrix array for fast access (direct) in other functions
        $rawArray: function() {
            return this.mat;
        },
        
        toArray: function() {
            var A = this.mat,
                arr = new Array(16);
            
            for (var i = 0; i< 16; i++) {
                arr[i] = A[i];
            }
            return arr;
        },
        
        log: function() {
            var A = this.mat;

            return [
                '[ ',
                    [A[0], ', ', A[4], ', ', A[8],  ', ', A[12], ', '].join(''),
                    [A[1], ', ', A[5], ', ', A[9],  ', ', A[13], ', '].join(''),
                    [A[2], ', ', A[6], ', ', A[10], ', ', A[14], ', '].join(''),
                    [A[3], ', ', A[7], ', ', A[11], ', ', A[15]].join(''),
                ' ]'
            ].join('\n');
        }

    } ); // module.Mat44 = Factory.obj()

    // friend methods ---------------------------------------------
    //

    var friends = {
        // multiply m by scalar and return m
        multipyByScalar: function(m, f) {
            var B = m.mat;
                 
            B[0] *= f; B[4] *= f; B[8]  *= f; B[12] *= f;
            B[1] *= f; B[5] *= f; B[9]  *= f; B[13] *= f;
            B[2] *= f; B[6] *= f; B[10] *= f; B[14] *= f;
            B[3] *= f; B[7] *= f; B[11] *= f; B[15] *= f;
            
            return m;
        },
        
        // transpose m and return it
        transpose: function(m) {
            var tmp,
                A = m.mat;
            
            tmp = A[4];
            A[4] = A[1];
            A[1] = tmp;
            
            tmp = A[8];
            A[8] = A[2];
            A[2] = tmp;
            
            tmp = A[9];
            A[9] = A[6];
            A[6] = tmp;
            
            tmp = A[12];
            A[12] = A[3];
            A[3] = tmp;
            
            tmp = A[13];
            A[13] = A[7];
            A[7] = tmp;
            
            tmp = A[14];
            A[14] = A[11];
            A[11] = tmp;

            return m;
        },
        
        // Row vector - Matrix multiply  ie. c = b[A]
        multiplyRow: function(rowVec, m) {
            var A = m.mat,
                b = rowVec.toArray();
            
            return module.Vec3(
                b[0]*A[0]  + b[1]*A[1]  + b[2]*A[2]  + b[3]*A[3],
                b[0]*A[4]  + b[1]*A[5]  + b[2]*A[6]  + b[3]*A[7],
                b[0]*A[8]  + b[1]*A[9]  + b[2]*A[10] + b[3]*A[11],
                b[0]*A[12] + b[1]*A[13] + b[2]*A[14] + b[3]*A[15]
            );
        },
        
        AffineInverse: function(m) {
            return Mat44(m).affineInverse();
        },

        affineInverse: function(m) {
            return m.affineInverse();
        },
        
        // negate m and return it
        negation: function(m) {
            var A = m.mat;
            
            for (var i =0; i< 16; i++) {
                A[i] = -A[i];
            }
            return m;
        },
                
        copyInto: function(from,to) {
            for (var i =0; i< 16; i++) {
                to.mat[i] = from.mat[i];
            }
            
            return to;
        },
        
        toArray: function(M, arr) {
            var A = M.mat;
            
            for (var i = 0; i < 16; i++) {
                arr[i] = A[i];
            }
            return arr;
        },

        $reader: function(m) {
            return function(iRow, jCol) {
                return m.mat[iRow + 4*jCol];
            }
        },

        // return the raw matrix array for fast access (direct) in other functions
        $rawArray: function(m) {
            return m.mat;
        }
        
    }; // friends{}

    for (var f in friends) {
        friends.hasOwnProperty(f) && (module.Mat33[f] = friends[f]);
    };
    
    var Mat44 = module.Mat44;


} )(mLIB);